A BEAM MODEL FOR HYDRAULIC FRACTURING 
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O . Abstract 



We investigate numerically the shape of cracks obtained in hydraulic fractur- 



ing at constant pressure using a square lattice beam model with disorder. We 
consider the case in which only beams under tension can break, and discuss the 
conditions under which the resulting cracks may develop fractal patterns. We also 
determine the opening volume of the crack and the elastic stress field in the bulk, 
quantities which are accessible experimentally. 
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Hydraulic fracturing is a technique of great importance in soil mechanics and 
is used systematically to improve oil recovery. In particular, it has a considerable 
importance for the functioning of geothermal wells An incompressible fluid, in 
general water, is pushed with pressure deep inside a solid, in the case of soil by 
injecting it into a deep perforation. The fluid penetrates into the solid by opening 
long cracks that radially go from the injection hole into the material. 

In two dimensions controlled experiments have been performed recently^. 
Water or air is pushed from above into the center of a circular Hele Shaw cell filled 
with clay. While on long time scales clay behaves like a fluid for high injection 
pressures it fractures like a solid (viscoelastic fluid). In this fracturing regime the 
resulting cracks display a disordered ramified structure which appears to obey self- 
similarity with a fractal dimension of 1.4 - 1.5. Compared to the fractal structures 
observed in Laplacian systems (dielectric breakdown, viscous fingering or diffusion- 
limited aggregation (DLA) ) not only the fractal dimension is lower but also the 
angles between branches are about three times larger (i.e. about 90°) and no tip 
splitting is observed. 

These interesting observations remain largely unexplained. While for the 
Laplacian case the underlying instabilities (Saffman- Taylor, tip splitting, side 
branching) have been heavily investigated! 3 ] and numerical clusters of tens of mil- 
lions of particles have been generated W, in the case of fracture much work still 
needs to be done. A recent stability analysis of the shape of a circular hole with 
either internal pressure or in a stretched membrane has shown that contrary to the 
common belief a large difference can be expected in the shape of a crack between 
the two cases^ 5 ]. The origin of this is the non-linear dependence of the growth 
velocity of the crack surface due to the threshold in cohesion force that must be 
overcome to break the material. In the Laplacian case these differences do not 
exist in the limit of large clusters^. 
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Numerically the breaking of a material from a central hole was first investi- 
gated by Louis and Guinea^ using a triangular network of springs and stretching 
the network radially on the outer boundary into the six directions of a hexagon. 
They observed fractal cracks having a fractal dimension close to that of DLA. 
These findings were studied subsequently by various groups^. It became clear 
that the patterns depend very much on the type of applied displacements (shear, 
uni-axial, radial). Since networks of springs need more extended structures to 
transfer momenta it turned out to be more efficient to use a beam model ^ which 
is no longer a central force model. For more details on the work that has been 
done we refer to the book of Ref. 10. 

While in previous works imposed displacements were applied only at the ex- 
ternal surface of the solid (Dirichlets boundary value problem), in this paper we 
introduce and study a model in which the imposed load represents a pressure that 
acts along the whole (inner) surface of the crack in a direction perpendicular to 
the surface (von Neumann boundary value problem). In this way, the point of 
application of the imposed load varies during the growth of the crack, a situation 
that more realistically describes the case of hydraulic fracturing. 

We consider the beam model (as defined in p. 232 of Ref. 10) on a square 
lattice of linear size L. Each of the lattice sites i carries three real variables: the 
two translational displacements Xi and yi and a rotational displacement <pi. Next 
neighbouring sites are rigidly connected by elastic beams of length l^. When a 
site is rotated ipi ^ then the beams must bend accordingly (Fig. la) and continue 
to tangentially form 90° angles with each other. In this way local momenta are 
taken into account. 

Since most materials typically crack under tension and in much less degree 
under compression we will assume that only beams that are under tension are 
allowed to break, i.e. to be irreversibly removed from the forces and momenta 
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balance equations^ 12 !. This means that the cohesion force against compression is 
actually infinite, i.e. the beams can be compressed but not broken by compression. 

At the place into which the incompressible fluid is supposed to be injected 
(center of the lattice) one vertical beam connecting sites i and j is removed. Since 
we want to simulate the loading of a crack by the injected fluid an invariant double 
force Fi and Fj conserving momentum is applied at the sites i and j pointing from 
the hole (removed beam) into the elastic bulk (Fig. lb). Thus, force densities are 

— * 

replaced by discrete force vectors F = (F x , F y ), 

F l = F (0,1), Fj = -F (0,1). (1) 

The forces for a broken horizontal bond are defined correspondingly. The acting 
pressure P is just the force F per beam thickness d and beam length / and its 
value is kept fixed during the fracture process. 

Each time a beam is removed a new force dipole is applied at its corresponding 
nearest- neighbor sites i and j which destroys the balance of forces existing previ- 
ously. Consequently the lattice must be relaxed and the new static equilibrium 
must be obtained again. This is done in our case using conjugate gradient' 13 ^. 
After the unique displacement field Ui corresponding to the new boundary condi- 
tions is determined one can decide which is the next beam to be broken. Due to 
this procedure we only consider the case of an instantaneous relaxation, i.e. we 
assume that the physical process of stress redistribution is much faster than the 
actual fracture process. 

For simplicity, only beams along the surface of the inner hole are considered in 
the breaking process^ 14 !. For each of these beams the force fy acting along its axis 
is determined and only if this fy is positive, i.e. a tension, the beam can be broken 
(active beams). In addition, since we are interested in random structures we need 
to incorporate some stochastic mechanism or disorder in the model. The disorder 
may be either quenched, i.e. the beams may have different thresholds for breaking, 
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or it may be annealed, i.e. one selects the next broken beam randomly according 
to some probability. In the following we will employ the second procedure. 

Among the active beams one is chosen with a probability pij , which is deter- 
mined in analogy to dielectric breakdown^ 15 ! to be 

Pij oc fl (2) 

if fij is larger than a certain cohesion force threshold / co h > and zero otherwise' 16 ' . 
Here, < r\ < oo is a parameter which determines the local breaking properties of 
the material. If rj — > oo only the most stretched beam breaks, i.e. the one with 
the largest fij. The opposite case corresponds to r\ = where all active beams 
can break with the same probability, i.e. to a situation in which heterogeneities 
are more important than tension forces. 

Let us start discussing our results for the case rj = 1 (similar results are 
obtained for r\ > 1) and / co h = 0. For the Laplacian case, the value r\ = 1 
corresponds to DLA with a fractal dimension df = 1.7 in two dimensions' 15 '. For 
the present beam model, however, topologically linear cracks develop (df = 1) as 
shown in Fig. 2a. This indicates that the distribution of is much narrower 
than for DLA and essentially only the most stretched beams (located at the tips 
of the crack) break. Thus, already for r\ = 1, only a single branch persist on large 
length scales. The colors in Fig. 2 represent the stress field. They show that the 
strongest gradients are indeed around the tips and that regions other than the tips 
are under constant compression of magnitude P. 

Also surprising is the result obtained when t] = 0, i.e. the case in which all 
beams under tension can be selected for breaking with the same probability. For 
the Laplacian case, one obtains compact clusters of spherical shape and fluctu- 
ations occur only at their surface! 15 '. In the present case, the resulting cracks 
(shown in Fig. 2b) display ramifications which persist up to scales comparable 
to the system size L, and it seems plausible that they might be fractal. The col- 
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ors show that the strongest gradients occur again around the crack tips and that 
the regions between crack branches are virtually under constant compression of 
magnitude P. One clearly sees the opening of the crack being larger on certain 
main cracks than on the side branches. Although only weakly visible the system 
of Fig. 2b is no perfect square anymore but has small bulges' 17 ]. 

To characterize quantitatively the shape of the cracks obtained when 77 = 0, 
we have plotted the number N of broken beams against the radius of gyration Rq 
of the crack' 18 !. The results are shown in the log-log plot of Fig. 3 for the case 
of free external boundaries. We see a powerlaw behaviour N ~ Rq over nearly 
two decades and the slope gives a fractal dimension of df = 1.56 ± 0.05. Similar 
results were obtained for the case of periodic boundary conditions. 

Our model is different with respect to the Laplacian one in that all the beams 
under compression cannot break. The fact that in our case the cracks become 
fractal when r\ = shows that only at very few points, namely at the crack tips, 
the crack surface is under tension and everywhere else under compression. 

It is also interesting to consider a fractal structure through its intrinsic metric 
by using as geodesies only the shortest paths that are entirely on the cluster. In 
this topological or "chemical" space ^ , one can calculate the number of broken 
bonds Ng that are found within a chemical distance £ ' 20 L We have obtained a 
powerlaw behavior Ng ~ £ de with dg = 1.4 ± 0.05, for the same cracks considered 
in Fig. 3, somewhat smaller than df. The interest in the "chemical dimension" 
dg is because for loopless structures, such as the cracks described here, both df 
and dg completely determine the scaling behavior of diffusion controlled transport 
phenomena I 21 l such as e.g. the diffusion of chemical substances which can be 
present in the fluid contained within the cracks. 

For a finite cohesion threshold / C oh > and r] = the cracks tend to grow 
more anisotropically since, initially, only the most stressed beams can break' 16 !. 
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However, much larger system sizes than the ones considered here are required in 
order to reach the asymptotic shape of the cracks when / co h > 0. This is because 
a finite ratio f co h/P introduces an additional length scale in the system I 22 !, and 
the scaling behavior obtained when / co h = is recovered only asymptotically. 

Our simulations are performed at constant pressure and the volume V of 
the crack opening which corresponds to the amount of fluid that has penetrated 
into the soil is a measurable variable. We measure this volume by taking into 
account the actual displacements of the lattice sites as the crack grows. These 
displacements are given from the elastic solution U{ = (xi, yi, cpi), and the volume 
elements AV^ connecting sites % and j are obtained simply as AV^ = (xi — Xj)ld 
for horizontal broken beams, with / and d being the beam length and thickness, 
respectively. Similar expressions hold for vertical broken beams. The total crack 
opening volume V is just the sum of all the volume elements AV^. 

In Fig. 4 we show the crack volume V for free and periodic boundary condi- 
tions as a function of the radius of gyration Rq in a log- log plot. The slope gives 
an exponent consistent with two, i.e. the spatial dimension. This agrees with the 
observation that a finite amount of fluid actually enters into the hole. When the 
crack approaches the boundary too much the curve has an artificial increase in 
slope due to boundary effects. For comparison, we also measured the volume of 
a deterministic straight crack as shown in Fig. 4. For such a crack in an infinite 
plate it is well known from elasticity theory that the crack opening volume is pro- 
portional to the squared crack length (in two dimensions) t 23 '. This shows that 
(1) although the crack surface (number of broken beams) is fractal the volume of 
the crack is not, (2) the relationship V ~ R? G is independent of the crack structure 
itself and (3) particular external boundary conditions play a minor role as long as 
the "typical length" Rq describing the crack structure is much smaller than the 
system size. 
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We have presented a discrete fracture model in two dimensions adapted to 
the conditions of hydraulic cracking of soils as typically used for instance to cre- 
ate additional heat exchange surface for geothermal energy recovery. Only tensile 
forces break the solid and the heterogeneities are considered to be dominant t 24 l 
Under these circumstances the crack patterns are fractal and we determined the 
fractal dimension to be less than that of DLA. Our results for r\ = reproduce 
some of the conspicuous features observed in the 2d experiments of Orleans for 
high injection rates and rigid pastes' 2 ': The crack patterns are more kinky than 
DLA and the fractal dimension is lower, their (large) error actually overlap. Thus, 
the basic assumptions of our model, namely that the material does not open un- 
der compression and that for 77 = the heterogeneities are more relevant than 
the actual value of the tension force, may describe the case of pastes with large 
clay/water ratios used in Orleans. 

We also have clear numerical evidence that at constant pressure the opening 
volume of the crack grows proportionaly to the area spanned by the crack, i.e. the 
square of its radius of gyration, independent of the crack structure itself. 

Our model should be considered as yet quite preliminary concerning the un- 
derstanding of real industrial hydrostatic fracturing. On one hand real soils are 
three dimensional and in fact we do not know whether in three dimensions the 
patterns are finger-like as in two dimensions or rather consisting of randomly at- 
tached fracture planes. The later case seems more likely has yet never been seen 
numerically up to now. Also the role of the heterogeneities must be investigated 
closer. In soils the disorder is quenched, i.e. the randomness (in breaking thresh- 
old, modulus, etc) are assigned before cracking starts. We have used "annealed" 
randomness: random numbers are used during the breaking process in order to 
select beams with a certain probability. Also realistic, e.g. Weibull distribu- 
tions, should be considered. Other important ingredients that should be taken 

8 



into account when trying to make the model more realistic are plasticity, pressure 
gradients and hydrodynamic effects in the fluid, stress corrosion and short time 
effects, like shock waves. 
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Figure Captions 

Fig. 1 The beam model on the square lattice, (a) A beam connecting sites i and j is 
shown to display the rotational displacements (p at both sites, (b) Diagram of 
the pair of external forces F (arrows) applied at sites i and j corresponding 
to the initially removed beam. 

Fig. 2a A typical crack obtained with a beam model for hydraulic fracturing at con- 
stant pressure when 77 = 1, on a square lattice of 200 x 200 sites and free 
external boundary conditions. The crack consists of 310 broken beams, and 
only beams under tension are allowed to break. The different colors represent 
the intensity of the hydrostatic stress field \fij\. The full range for the stress 
field is linearly mapped onto twenty color-circle cycles starting at Magenta 
going through Blue, Cyan, Green, Yellow and ending at Red. 

Fig. 2b Same as in Fig. 2a for 77 = 0, on a square lattice of 150 x 150 sites. This crack 
consists of 680 broken beams. 

Fig. 3 Plot of the number of broken beams N as a function of the radius of gyration 
Rq of the cracks obtained when 77 = and free boundary conditions. Averages 
over 27 samples were performed. The inset shows the successive slopes df of 
the data, and indicate an average value df = 1.56 ± 0.05. 

Fig. 4 Plot of the crack opening volume V as a function of the radius of gyration Rq, 
obtained by averaging over 12 samples using free boundary conditions (trian- 
gles) and averaging 4 samples using periodic boundary conditions (squares). 
For comparison, we also show the results for a straight crack (circles) obtained 
in the case 77 — > 00 under free boundary conditions. In all three cases we find 
the behavior V ~ R 2 G . 
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